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Abstract. We argue that the main mechanism for condensate collapse in an 
attractive Bose-Einstein condensate is the loss of coherence between atoms a 
finite distance apart, rather than the growth of the occupation number in non 
condensate modes. Since the former mechanism is faster than the latter by a 
factor of approximately 3/2, this helps to dispel the apparent failure of field 
theoretical models in predicting the collapse time of the condensate. 
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1. Introduction 

The so-called Bose Nova experiment on the collapse of a Bose-Einstein condensate 
with attractive interactions [U El [3] has opened up a fascinating window in the far 
out of equilibrium behavior of these systems. The experiment has been analyzed 
from a number of perspectives [H [7j [9l [10] and it is fair to say that we have 
a good qualitative understanding of the phenomenon. However, at the quantitative 
level certain anomalies persist. 

In this paper we shall deal with the apparent failure of existing models in 
predicting the collapse time scale t c for the condensate, in the regime where the 
scattering length a is just below the critical value — a c . In [7J the scaling law 



l£l _ j 



-1/2 

(1) 



is proposed, which fits well the experimental results. However, the proportionality 
constant is not derived. The authors of [Jj claimed that the proper proportionality 
constant could be derived from a complete field theoretic calculation, but when the 
calculation was actually done [HI H21 H3j , it failed to produce a satisfactory prediction. 

In this paper we shall present a qualitative analysis of the collapse time for 
a condensate trapped in a flat box [5] with periodic boundary conditions. Unlike 
previous analysis, we shall assume that the total number of particles in the condensate 
remains fixed [131 [TSl HH UJ HH] • Under these constraint, the condensate occupation 
number is properly defined as the greatest eigenvalue of the one particle density matrix 
(to be defined below) [19) . Given the assumed geometry, the corresponding eigcnmodc 
is necessarily homogeneous, so the eigenvalue is just the integral of the one particle 
density matrix with one argument fixed, and the other ranging over the confining box. 

Since the overall normalization of the one particle density matrix is determined by 
the overall density of the gas (see below), the most important factor in the evolution of 
the condensate occupation number is how fast the density matrix falls off, measuring 
the degree of coherence among atoms at finite distances. We shall argue below that the 
one particle density matrix is approximately Gaussian with a variance which decays 
in time as exp{— jt}, with 




(2) 



where 



r 1 - ^ (3) 
T " ML 2 [6) 

is the frequency of the first excited states for a particle in the box; here M is the 
mass of an atom and L is the size of the box. Therefore, after integrating over the 
three dimensional box we obtain that the condensate occupation number decays as 
exp {— 3jt}. 

The expectation number in the first excited state, as computed from the Gross 
- Pitaievskii equation, the Hartree - Fock - Bogoliubov or the Popov approximation 
would grow only as exp{27i}. Therefore, condensate collapse from loss of coherence 
between atoms is faster than the estimate from loss of particles to excited modes by 
a factor 3/2. For comparison, note that a detailed calculation of the collapse time for 
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a = — 10a c yields a predicted value of 10ms against an experimental value of (6 ± 1) ms 
[13] . Therefore a factor of three halves goes a long way to solve the existing puzzle. 

This paper is organized as follows. In next section we present the model of a cold 
trapped Bose gas, introduce a suitable set of density and phase variables and solve 
the Heisenberg equations in the linearized approximation. In the following section, 
we apply these results to derive the evolution of the condensate particle number and 
thereby our main result. In Section 4 we compare this results to the particle number 
conserving, Hartree-Fock Bogoliubov and Popov approaches. We close the paper with 
some brief final remarks, and give supplementary technical details in the Appendix. 



2. The model 



The idea is to analyze the Bose Nova experiment with the tools we have developed 
to handle the Mott transition in [20] . The starting point is a second - quantized field 
operator ^> (x, t) which removes an atom at the location x at times t. It obeys the 
canonical commutation relations 



*(x,t),*(y,t) 



= 



The dynamics of this field is given by the Heisenberg equations of motion 
d 



where H is the Hamiltonian. 
the field operator 



(4) 
(5) 



(6) 



The theory is invariant under a global phase change of 



e~ lf> ¥ (7) 

The constant of motion associated with this invariance through Noether theorem is 
the total particle number. 

To progress further, we need a specific model for the atom-atom interactions. 
In principle, we should specify the atom-atom interaction potential. However, in 
many applications it is enough to know the cross section a for low energy spherically 
symmetric scattering of two identical bosons. We introduce the scattering length a 
through a = 87ra 2 , where the factor 8ir involves both integration over all scattering 
angles and Bose enhancement factors. We shall adopt as model atom-atom interaction 
a contact potential U5 (x) . This is expected to be a good approximation as long as 
the distance between atoms is much greater than both the scattering length and the 
distance out to which the fundamental atom-atom interaction is important |21j . To 
reproduce the right scattering length we need U — 4Trh 2 a/M, where M is the mass of 
the atoms. A positive value of a means a repulsive interaction; we adopt the convention 
that an attractive interaction is described by a negative value of a. 

Assuming a contact atom-atom potential we get the Hamiltonian 



H = < * T #* 



U 



The single-particle Hamiltonian H is given by 

h 2 



H ^ = "^7 V * + V t™p (X) * 
Z1V1 



(8) 



(9) 
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V trap denotes a confining trap potential. Then the Heisenberg equation of motion 

d » ~ - » . 
i%— * = + C/*^ 2 (10) 
at 

is also the classical equation of motion derived from the action 

S = ( dtdxim*^ - J dtn (11) 

placing hats everywhere. For simplicity we shall replace the trap potential by a flat 
bounding box of volume V — L 3 with periodic boundary conditions. Yurovsky has 
demonstrated that this is enough for a qualitative treatment of the Bose Nova [5] . We 
also assume that we have a finite total number of particles TV, which remains fixed 
through the evolution (that is, there is no particle loss to the environment). 

2.1. Density and phase variables in the CTP formulation 

To analyze further this model we shall adopt density-phase variables [22, 23 . These 
variables have been extensively used to study dynamical problems, including the Mott 
transition [24 . This will set the stage for a further canonical transformation to a more 
convenient set of degrees of freedom, to be carried out in the next Section. 

In the path integral representation, quantum amplitudes are given in terms of 
functional integrals over complex fields \& and ^ associated to the destruction and 
creation operators. Our starting point is the Madelung representation [22l [23] 



*(x,t) = [exp-*p(x,i)] VplxJ) (12) 
tft (x,t) = Vp(*,t) [exp^(x,i)]. (13) 

In the canonical formalism, the fields p and ip become operators with commutation 
relations 

[p(x,t),0(y,t)] = -tf(x-y), (14) 

Within the path integral we allow the phases ip to take all real values, and therefore 
so do the conjugated density operators p [501 US]- This makes the square roots in ([12"]) 
and ([13]) problematic. It is best to adopt a new set of variables where square roots do 
not appear, as we shall do presently. For further discussion of density-phase variables 
in continuum theories see [26] , 

We adopt the formalism developed in [20] to describe the transition from the 
superfluid to the Mott insulator state in an optical lattice. To compute expectation 
values, we shall use the closed time-path formalism, where we choose the independent 
variables as follows. In the first branch, we define a new (complex) variable x 1 ( x > t) 
from 

^(x,*) =exp[- lX 1 (x,i)] (15) 
* 1 t(x,t)-p 1 (x,i) expfix^x,*)] (16) 
On the second branch we write instead 

* 2 t(x,t)=exp[ J x 2t (x,t)] (17) 
* 2 (x, t) = exp Hx 2t (x, t)] p 2 (x, t) (18) 
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In the canonical formulation, the fields x an d P become operators with commutation 
relations [5U] 

[p(x,t),*(y,t)]=-W(x-y), (19) 
The dynamics of these operators is given by the Hamiltonian 

H(p, X ) = ydx|A_(pVx"*Vp)Vx+^/5(/3-l)| (20) 

plus the necessary terms to enforce a fixed total particle number [20 . Observe that 
in the new variables, the action is explicitly analytical. 

We now split all variables into a homogeneous and an inhomogeneous part. 

p(x,t) = n + f(x,f) (21) 
r(x,*) = ^r p (t)/ p (x) (22) 

where the / p are plane waves 

fp ( x ) = yrj2 exp {*P x / n } ( 23 ) 

and the allowed values of the components p^, /i = 1 — 3, of the momentum p are 
integer multiples of 2~k%/ L, and similarly 

*(x,*) = ^ + X;*p(*)/p( x ) (24) 

Observe that the homogeneous part of the density operator is constrained to be the 
c-number n = N/V, and the homogeneous part of the phase is a collective coordinate 
[27] which couples only to the homogeneous density. It does not affect the dynamics 
of the inhomogeneous modes. 

Consider the lowest order theory which is obtained by keeping only the "free" 
quadratic part of the Hamiltonian 



where v p = p 2 /M , p = |p|. The Heisenberg equations of motion are 



(25) 



d • —iv„ ~ 

k di Xp = ~J^ Xp + Ufv (26) 
- h J t f P = nv v^v + — Y- f P ( 27 ) 



Where 



4Trh 2 Na 

Un =^njr ( 28 ) 

In the Bose Nova scenario, we have U = if t < 0. Therefore the frequencies are 
fiLOp = Vp/2. If we call A p the destruction operator which kills the initial state, then 

r p (0-) = (-i) Up - Al p ] (29) 
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x P (o-) 



-.A r 



(30) 



For t > 0, we have U < instead, and 



X p (t) = -L | cos [w p t] - * (y + f 7 ") 



it/^i p 



with the dispersion relation 



1 



sin [wpt] 
huj r . 



(31) 



(32) 



3. The one-body density matrix 

We may now turn to computing the one-body density matrix 

a (x, y, t) = (*t ( X; t ) $ ( y , t )) = ( exp < [^* (X) t) _ x i (yj t) ] j (33) 

In the last term, the 1,2 supcrindcx indicates closed-time-path ordering: operators 
with a 2 superindex always go to the left of operators with a 1 superindex. Observe 
that in our variables, the observable to be computed is a pure exponential: there 
are no square roots to be developed. This is the whole point of introducing the new 
variables. 

As in the previous Section, we separate the variables \ 2 * an d X 1 m to their 
homogeneous and inhomogeneous parts. Observe that the homogeneous terms may 
affect the overall normalization of the one particle density functional but not its shape. 
The overall normalization, on the other hand, is determined by the requirement that 
a (x, x, t) = n. So we may simply continue to disregard the homogeneous terms. 

Since we have restricted ourselves to a Hamiltonian which is quadratic in the 
inhomogeneous modes, we may use the Wick theorem result 



(^) = (1>Exp{y<^)} 



(34) 



with 



A = A(x,y)= X 2 *(x,i)-x 1 (y,<) (35) 
Decomposing the field operators in modes, with due attention to the closed time path 
ordering, we obtain 

1 



(A 2 ) = const. + 2^2 
Using the decomposition (|3Tj) 



V 



f-v (x) / P (y) 



(A 2 ) = const. + 



4(Z7n) 



p#0 



P (x - y) 
2h 



2 r 



sin [u) p t] 



1 2 



hu! r , 



(36) 



(37) 



To continue, we consider only the contribution from the unstable modes. The condition 
for instability is U < with \Un\ > v p /A. Since the lowest possible nontrivial value 
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of p is h/L, we get the critical scattering length as a c = ttL/AN. For a close enough 
to the critical value, the six modes with L 2 p 2 = h 2 are the only unstable ones. Their 
frequency is u> — —17, where 7 is given in ^ above. Setting y = and x = |x|, we 
get 



E 

p=h/L 



I sir 



rpx- 


r 


sin [ujpt] 


2 

- 2 


7r sinh [yt] 


"(z) 


i2h- 




Tiujp 




fvy 





therefore 



a (x, t) — n Exp ■ 



2%U n sinh [jt] 



N l / 2 hjL 



(38) 



(39) 



The condensate occupation number N c is obtained by integrating over x, so, once the 
Gaussian approximation becomes valid 



N r cx e 



-37i 



(40) 



We therefore obtain the same scaling law as in [7], but the coefficient is 3/2 times 
larger. As noted in the Introduction, this correction is enough to account for the 
anomaly observed in [T3] . 



4. Comparison to other approaches 

In this Section we will compare the result above for the one-particle density matrix 
with other approaches in the literature, namely the particle-number conserving (PNC) 
formalism and the Hartree-Fock-Bogoliubov (HFB) and Popov approximation. We 
shall not discuss the so-called truncated Wigner approximation, but refer the reader 
to the detailed treatment in [13]. See [28] for further details on these approaches. 



4--1- The equations of motion in the PNC approach 

The PNC formalism [TJl HH EH El HE] is usually presented as an expansion in inverse 
powers of the total particle number N. In preparation for this, it is convenient to scale 
the interaction term, writing U = u/N. 

The basic insight of the PNC approach is that if the total particle number 
remains constant, then each particle above the condensate corresponds to a hole in 
the condensate, so we may speak of particle-hole (PH) pairs. 

Let us consider the expansion of the field operator in plane waves 

*(x,i) = 5> p (i)/ p (x) (41) 
p 

ao reduces the number of particles in the condensate by one. Following Arnowitt and 
Girardeau, let us introduce the operator 

(3 = J- «o = a Q ^= (42) 

VAq + i Vn 

where 

N = N-Y, a P «P (43) 

p^O 
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is the condensate number Heisenberg operator. Observe that for a number eigenstate 
f3\N ) — |JV"o — 1) unless Nq = 0, in which case (3 |0) = 0. Therefore (3 preserves 
the norm for all states orthogonal to the state with no particles in the zeroth mode 
(which is much stronger than not having a condensate). If there is a condensate, any 
physically meaningful state will satisfy this requirement, and (3 may be considered a 
unitary operator, with inverse 

/J t = ^aS = 4—L= (44) 

We now introduce the destruction operator of a PH with the particle in mode p 

A P - /3V (45) 

If we consider the /3's as unitary, then the A's satisfy bosonic canonical commutation 
relations. This relationship may be inverted 

a P = (3X P (46) 

The number of particles in a given mode is equal to the number of PH 

apa p = A p Ap (47) 

We write the field operator restricted to the subspace with a well defined total 
number of particles N as ^ = y/~N (3<j) 

4> = 0o + ^A (x, t) - [Sn (t)} 0o (48) 

where for a homogeneous condensate we must have 4>q = V^ 1 ^ 2 

A(x,i) = ^Ap(t)/ p (x) (49) 

5n(t) = J d 3 x A f A (50) 



F (x) = 2N 



1 - 











0(N-') (51) 



Within our approximations (3 commutes with <j>. To lowest order in AT 1 , A evolves 
according to 

= -ihX, t +H\ + Un(X + A^ + O (a^ 1/2 ) (52) 

(see Appendix) 
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4-2. The HFB and Popov approximations 

Before proceeding to compute the one particle density matrix in the PNC approach, let 
us show that the HFB and Popov approximations give essentially equivalent results. 

The HFB and Popov approximations are implementations of the symmetry 
breaking approach to condensation, where the formation of a BEC is associated to 
the spontaneous breaking of the J7 (1) symmetry (0 [29]. The field operator develops 
a c-number expectation value, which by traslation symmetry may depend only on time 
(*) = e -«W$ (t) (53) 

More generally 

* = e- l0 W [$ (t) + ip] (54) 

In the HFB approach, we use this decomposition in the Heisenberg equations 
of motion, where we also replace products of two fluctuation operators by their 
expectation value, and use the so-called Hartree approximation. 

ip^ip 2 ~ 2hip + rhipi (55) 

where 

n = (56) 

m = (if 2 ) (57) 

The Heisenberg equations decompose into equations for the mean fields and equations 
for the fluctuations 

ih^-® + -q<5> = ?7$ 3 + 2Un<S> + Urh>5> (58) 



dt 
d 

i^gli> + # = Hi]) + 2U ($ 2 + n) ip + U(<P 2 + m) ifi (59) 



where 

,.«* (60) 

The HFB approximation has the serious drawback that it is not gapless, and therefore 
it is hardly reliable in a problem such as the Bose Nova, which depends critically on 
the behavior of long wavelength modes. The Popov approximation overcomes this 
problem by further neglecting rh. Then we obtain 

r? = U<S> 2 + 2Un (61) 

and the fluctuation equation becomes 

g 

ih—Tf> = H4> + U<5> 2 (iJj + ^) (62) 
at x ' 

Under this approximation the number of particles in the condensate remains constant. 
This may be avoided by including explicitly the effect of particle loss through three 
body recombination. However, the final results are robust against these terms 
[H QUI [13], and we shall not consider them in detail. On the other hand, the total 
number of particles is not conserved. 

If we assume that the temperature is effectively absolute zero, then <£> 2 = n 
initially and remains close to it until much later in the collapse; the effect of finite 
temperature is discussed in [13j and it is seen to be minor. If we just replace 
<J> 2 = n, the Popov equation for the fluctuations reduces to the PNC equation for 
the inhomogeneous modes (|52p . This approximation gives a reasonable description of 
early jet and burst formation [7 , so it may be considered reliable. 
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4-3. The one particle density matrix in the PNC approach 

We now return to the calculation of the one particle reduced density matrix. 



a(x,y,i) = (¥( X ,t)^(y,t) 



= n |l__[^„)_y( A t ( x ,t)A(y,t)) 
Decomposing in modes, we get 

<7(x,y,i) = rJ 1- 



p#0 



- /-p (x) / P (y) 



(A P Ap) 



Each mode evolves according to 



(63) 



(64) 



(65) 



dt 2 ' p 

The dispersion relation is given by (|32|) . At t = 0, A must destroy the physical state, 
so A p (0) = e l<Pp A p for some phase <^ p . From the equation of motion we derive the 
initial velocity 

lh llt (0) = ^ e ' VpA P + Un ( cIVpA p + e _ ^ p Al p ) (66) 

Therefore 



A P (t) 



z?7 n 



cos [w p t] — * ("2" + ^ 



sin [w p rj] 



(67) 



This equation and (|31[) show that 

(A p A p )=n<X p X p ) (68) 

and therefore the PNC result is just the first term in the expansion of our earlier 
result ([36]) in inverse powers of A 1 / 2 . 

Indeed, the representations of the field operators (fT5|) and are equivalent, 
to next to leading orden in TV -1 / 2 , provided we identify e l(Pp = —i and 
expj-iXo/yVaJ = „i/2/J. 

5. Final remarks 

After this point, it only remains to comment on the reasons why this proposal works. 

From the formal point of view, our expression for the one particle reduced density 
matrix is seen to agree with the perturbativc implementation of the particle number 
conserving approach to next to leading order. This agreement suggests that, more 
generally, our approach implements a resummation of the PNC expansion. A key 
feature is that the we use variables that keep the exponential structure of the one 
particle density matrix. Therefore, the method suggested amounts to a perturbative 
evaluation of the exponent, but it is non perturbative with respect to the final result. 
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This formal advantage of the proposed method correlates with a shift in the 
physical emphasis, from particle creation in the excited modes to the loss of coherence 
among atoms. Comparing this to other formal studies of decoherence, it comes as no 
surprise that the later process is faster than the former J3U1 131] • 

We submit this minor contribution with the expectation that it will help clear 
the way to a full quantitative understanding of this fascinating experiment. 
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Appendix: Derivation of (152[) 

The idea is to seek a solution of the Heisenberg equations of motion for 'J where (3 
and the A's have developments in inverse powers of N . Define a q-number chemical 
potential /t from 

= (69) 
' dt h K ' 



We have 



d 

ih—(j)={H-pL)(j> + u<iJ(t> 2 (70) 



We then find 



fl<po + U(j>Q 



[-ih\, t + {H - A) A + ucpl (2A + At)] 

V N 



+ 0(N~ 1 ) (71) 
Taking the expectation value we find 

= - (A) 0o + u<p 3 - — (p,X) 

V -/V 



+ 0(N^) (72) 
Recall that jj, is Hermitian. So we may decompose this equation into 

= - (A) 0o + u4>l - -4= <AA + At A) + O {N- 1 ) (73) 

and 

= (AA - At A) + O (N- 1 ) (74) 



Subtracting the expectation value from the Heisenberg equation, we get 
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o = m - a) <i> 



-L [-ihX, t + (if - £) A + u^l (2A + At)] + -L (£A) 



+ 0(N- 1 ) (75) 
and from l[7g]l .([7i |l and fl73|) we get 

fi = <fi) + (N- 1 ) ~^ = Un (76) 

Observe that this implies 

(fiX) = O (iV- 1 / 2 ) (77) 
The equation for A simplifies into (j5l2"|) 
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